PCA
library(readxl)
mydata <- read_xlsx("./Proteins.xlsx")
mydata <- as.data.frame(mydata)
head(mydata, 3)## Country Alcohol Cereals Eggs Fish Fruits Meat Pulses Spices
## 1 Albania 0.1840 14.2331 1.8069 0.6274 1.2757 7.6582 1.4555 0.0000
## 2 Austria 0.6150 13.1508 2.1456 1.8358 0.6242 13.6549 0.2405 0.1341
## 3 Belarus 0.3401 13.7248 2.2901 2.4548 0.4410 14.2402 0.0053 0.0213
## Treenuts Vegetables
## 1 0.2677 3.2456
## 2 0.3561 1.5259
## 3 0.3454 2.8587Description:
- Country
- The remaining variables represent the average amount of protein consumed by people in each country
mydata_PCA <- mydata[c(2:11)]
rownames(mydata_PCA) <- mydata$Country
desc_stats <- data.frame(
N = apply(mydata_PCA, 2, function(x) sum(!is.na(x))),
Missing = apply(mydata_PCA, 2, function(x) sum(is.na(x))),
Mean = colMeans(mydata_PCA, na.rm = TRUE),
Variance = apply(mydata_PCA, 2, var, na.rm = TRUE),
Min = apply(mydata_PCA, 2, min, na.rm = TRUE),
Max = apply(mydata_PCA, 2, max, na.rm = TRUE)
)
round(desc_stats, 2)## N Missing Mean Variance Min Max
## Alcohol 38 0 0.46 0.04 0.16 1.17
## Cereals 38 0 14.23 6.87 9.38 20.47
## Eggs 38 0 1.73 0.22 0.79 3.09
## Fish 38 0 3.03 3.36 0.63 9.40
## Fruits 38 0 0.68 0.04 0.39 1.28
## Meat 38 0 12.17 5.11 6.71 15.73
## Pulses 38 0 0.81 0.30 0.01 2.28
## Spices 38 0 0.14 0.05 0.00 1.35
## Treenuts 38 0 0.35 0.04 0.03 0.89
## Vegetables 38 0 1.99 0.60 0.98 4.37R <- cor(mydata_PCA)
library(GGally)
ggpairs(
mydata_PCA,
upper = list(continuous = wrap("points", size = 1)),
lower = list(continuous = wrap("cor", size = 2.5), colour = "black"),
diag = list(continuous = wrap("densityDiag", size = 0.6))
) +
theme(
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5),
strip.text.x = element_text(size = 6),
strip.text.y = element_text(size = 4),
axis.title = element_blank())## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the GGally package.
## Please report the issue at
## <https://github.com/ggobi/ggally/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning
## was generated.
# Convert correlation matrix to a distance matrix
dist_vars <- as.dist(1 - abs(R))
# Hierarchical clustering of variables
hc_vars <- hclust(dist_vars, method = "average")
# Reorder the correlation matrix according to the clustering
R_reordered <- R[hc_vars$order, hc_vars$order]
library(GGally)
ggpairs(R_reordered,
upper = list(continuous = wrap("points", size = 1)),
lower = list(continuous = wrap("cor", size = 2.5), colour = "black"),
diag = list(continuous = wrap("densityDiag", size = 0.6))) +
theme(axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5),
strip.text.x = element_text(size = 6),
strip.text.y = element_text(size = 5),
axis.title = element_blank())
library(psych)
cortest.bartlett(R, n = nrow(mydata_PCA)) ## $chisq
## [1] 107.6373
##
## $p.value
## [1] 4.766781e-07
##
## $df
## [1] 45det(R) ## [1] 0.03769258library(psych)
KMO(R) ## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R)
## Overall MSA = 0.6
## MSA for each item =
## Alcohol Cereals Eggs Fish Fruits Meat
## 0.41 0.55 0.44 0.46 0.50 0.80
## Pulses Spices Treenuts Vegetables
## 0.84 0.65 0.64 0.78library(FactoMineR)
components <- PCA(mydata_PCA,
scale.unit = TRUE,
graph = FALSE)
library(factoextra)
get_eigenvalue(components) ## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.1121916 31.121916 31.12192
## Dim.2 1.9695156 19.695156 50.81707
## Dim.3 1.1420515 11.420515 62.23759
## Dim.4 0.9863102 9.863102 72.10069
## Dim.5 0.6493707 6.493707 78.59440library(factoextra)
fviz_eig(components,
choice = "eigenvalue",
main = "Scree plot",
ylab = "Eigenvalue",
xlab = "Principal component",
addlabels = TRUE)## Warning: This FactoMineR PCA result contains only 5 eigenvalues and
## does not include the complete spectrum. Refit the PCA with a larger
## `ncp` before drawing a complete scree plot.
library(psych)
fa.parallel(
mydata_PCA,
fa = "pc",
cor = "cor",
n.iter = 1000,
sim = TRUE,
error.bars = TRUE
)
## Parallel analysis suggests that the number of factors = NA and the number of components = 2library(FactoMineR)
components <- PCA(mydata_PCA,
scale.unit = TRUE,
graph = FALSE,
ncp = 2)
components## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 38 individuals, described by 10 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"print(components$var$cor) ## Dim.1 Dim.2
## Alcohol -0.28418398 -0.62609438
## Cereals 0.70595796 -0.41718411
## Eggs -0.36437452 -0.34783780
## Fish -0.55347719 0.47844021
## Fruits 0.38926370 0.68325045
## Meat -0.74337993 0.03609410
## Pulses 0.57773132 0.36995101
## Spices 0.66082334 -0.24100811
## Treenuts -0.02573995 0.62422247
## Vegetables 0.78656544 -0.02909236print(components$var$contrib) ## Dim.1 Dim.2
## Alcohol 2.5949731 19.90307538
## Cereals 16.0136874 8.83682173
## Eggs 4.2660865 6.14319236
## Fish 9.8431277 11.62240247
## Fruits 4.8687949 23.70284231
## Meat 17.7564172 0.06614744
## Pulses 10.7247085 6.94910716
## Spices 14.0315105 2.94919779
## Treenuts 0.0212887 19.78424008
## Vegetables 19.8794054 0.04297328library(factoextra)
fviz_pca_var(components,
repel = TRUE)
library(factoextra)
fviz_pca_biplot(components,
repel = TRUE) 
library(factoextra)
fviz_pca_biplot(
components,
geom.ind = c("point", "text"),
geom.var = c("arrow", "text"),
label = "all",
repel = TRUE,
col.ind = "black",
col.var = "darkgreen",
pointsize = 1,
labelsize = 4,
alpha.ind = 1,
title = "PCA biplot") +
theme_minimal(base_size = 14)
library(factoextra)
library(ggrepel)
library(ggplot2)
ind <- as.data.frame(get_pca_ind(components)$coord)[, 1:2]
vr <- as.data.frame(get_pca_var(components)$coord)[, 1:2]
names(ind) <- c("x", "y"); ind$name <- rownames(ind)
names(vr) <- c("x", "y"); vr$name <- rownames(vr)
var_pct <- function(res) {
ev <- try(factoextra::get_eigenvalue(res), silent = TRUE)
if (!inherits(ev, "try-error") && !is.null(ev)) {
if (is.data.frame(ev) || is.matrix(ev)) {
j <- grep("variance.percent|percentage of variance|^variance",
colnames(ev), ignore.case = TRUE)
if (length(j)) return(as.numeric(ev[, j[1]]))
return(100 * as.numeric(ev[, 1]) / sum(as.numeric(ev[, 1])))
}
if (is.numeric(ev)) return(100 * ev / sum(ev)) # <- tvoj primer
}
if (!is.null(res$sdev)) return(100 * res$sdev^2 / sum(res$sdev^2)) # prcomp / princomp
if (!is.null(res$eig)) { e <- if (is.matrix(res$eig)) res$eig[, 1] else res$eig
return(100 * e / sum(e)) } # FactoMineR / ade4
if (!is.null(res$Vaccounted)) return(100 * res$Vaccounted["Proportion Var", ]) # psych
if (!is.null(res$values)) return(100 * res$values / sum(res$values))
rep(NA_real_, 2)
}
pv <- var_pct(components)
lab <- function(i) if (is.na(pv[i])) paste0("PC", i) else sprintf("PC%d (%.1f %%)", i, pv[i])
## --- 2. skupni faktor za pusice (ohrani kote med vektorji) -----
sc <- 0.85 * min(max(abs(ind$x)) / max(abs(vr$x)),
max(abs(ind$y)) / max(abs(vr$y)))
vr$xend <- vr$x * sc
vr$yend <- vr$y * sc
nud <- 0.10 # 10 % dolzine vektorja navzven
vr_layer <- rbind(
data.frame(x = vr$xend, y = vr$yend, name = vr$name,
nx = vr$xend * nud, ny = vr$yend * nud),
data.frame(x = ind$x, y = ind$y, name = "", nx = 0, ny = 0))
ind_layer <- rbind(
ind[, c("x", "y", "name")],
data.frame(x = vr$xend * (1 + nud), y = vr$yend * (1 + nud), name = ""))
ind_layer$nx <- 0; ind_layer$ny <- 0
pop_var <- data.frame(name = c("Fruits", "Pulses"),
dx = c( 0.35, 0.30),
dy = c( 0.30, -0.25))
i <- match(pop_var$name, vr_layer$name)
vr_layer$nx[i] <- vr_layer$nx[i] + pop_var$dx
vr_layer$ny[i] <- vr_layer$ny[i] + pop_var$dy
## --- 4. graf ---------------------------------------------------
p <- ggplot() +
geom_hline(yintercept = 0, colour = "grey80", linewidth = 0.3) +
geom_vline(xintercept = 0, colour = "grey80", linewidth = 0.3) +
geom_point(data = ind, aes(x, y), shape = 16, colour = "grey60", size = 1.1) +
# imena drzav: majhna, siva, z belo obrobo crk -> berljiva tudi cez mrezo
geom_text_repel(data = ind_layer, aes(x, y, label = name),
nudge_x = ind_layer$nx, nudge_y = ind_layer$ny,
size = 2.6, colour = "grey10",
segment.colour = "grey70", segment.size = 0.2,
min.segment.length = 0.2,
box.padding = 0.2, point.padding = 0.25,
bg.colour = "white", bg.r = 0.12,
max.overlaps = Inf, seed = 1) +
# spremenljivke: crne pusice
geom_segment(data = vr, aes(x = 0, y = 0, xend = xend, yend = yend),
arrow = arrow(length = unit(6, "pt"), type = "closed"),
colour = "black", linewidth = 0.5) +
# imena spremenljivk: manjsi krepki napis v belem okvircku, na konici pusice
geom_label_repel(data = vr_layer, aes(x, y, label = name),
nudge_x = vr_layer$nx, nudge_y = vr_layer$ny,
fill = "white", colour = "black", fontface = "bold",
size = 3.1, label.size = 0.25,
label.padding = unit(0.15, "lines"), label.r = unit(0.08, "lines"),
box.padding = 0.2, point.padding = 0,
segment.colour = "black", segment.size = 0.3,
min.segment.length = 0, # crta se narise vedno
force = 0.6, force_pull = 1.6, # oznaka ostane pri pusici
max.overlaps = Inf, seed = 3) +
coord_equal(clip = "off") +
scale_x_continuous(expand = expansion(mult = 0.08)) +
scale_y_continuous(expand = expansion(mult = 0.08)) +
labs(title = "PCA biplot",
x = lab(1), y = lab(2)) +
theme_minimal(base_size = 12) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "grey93", linewidth = 0.3),
plot.title = element_text(face = "bold"))
p
ggsave("pca_biplot.png", p, width = 8, height = 5, dpi = 600, bg = "white")